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FOREWORD 


Target  motion  analysis  (TMA)  is  a critical  function  in  naval  combat  con- 
trol systems.  In  many  situations,  bearing  is  the  only  information  available  and 
bearings-only  target  motion  analysis  is  the  foundation  of  the  TMA  process. 

This  report  is  the  result  of  research  directed  toward  the  estimation  of  target 
position  and  motion  parameters  via  bearing  measurements  generated  by  a single 
moving  observer.  However,  the  modeling  and  theory  presented  are  applicable 
to  the  general  case  of  multiple  moving  or  stationary  observers.  The  estima- 
tion algorithm  can  also  be  applied  to  navigational  systems  in  a situation  where 
a vehicle  observes  the  bearings  to  transmitters  of  known  location.  The  theory 
presented  in  this  report  can  be  extended  to  include  the  processing  of  multiple 
sensor  information  when  available  to  improve  the  performance  and  capability 
of  TMA  in  naval  applications.  Implementation  of  this  estimation  algorithm 
should  provide  significant  improvement  in  the  convergence  properties  of  auto- 
matic TMA  solutions. 
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POSITION  AND  VELOCITY  ESTIMATION 
VIA  BEARING  OBSERVATIONS 


INTRODUCTION 


The  problem  of  estimating  the  position  and  velocity  of  a body  traveling  on 
a straight-line  course  from  a history  of  bearing  measurements  is  considered  in 
this  report.  The  body  is  assumed  to  emit  an  appropriate  signal;  and,  ir  the 
general  two-dimensional  problem,  the  bearings  are  viewed  from  observers  con- 
fined to  a plane  that  includes  the  target.  The  situation  is  illustrated  in  figure  1 
wnere  it  is  desired  to  estimate  the  position  vector  rp  = (rpxrpy)'  velocity 
vector  Vp  = (VpxVpy)'  of  the  observed  body  from  noise-corrupted  measure- 
ments of  the  bearing  $.  The  inverse  situation  where  the  body  observes  the 
bearings  to  transmitters  of  known  location  represents  a system  capable  of  pro- 
viding navigational  information. 


y 


Figure  1.  Geometry  for  General  Two-Dimensional  Bearings-Only 
Motion  Analysis  Problem 

Defining  xp  = (rpxrTy  VpxVpy)'  as  the  target  state,  xQ  = (r0xrOyvOx 
VQy)’  as  the  observer  state,  and  xr  = xp  - xq  as  the  relative  target  state, 
the  problem  of  determining  xr  has  been  posed  as  a nonlinear  estimation 
problem.  A common  method  employed  to  obtain  a solution  has  been  the  extended 
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Kalman  filter,  although  this  method  has  experienced  convergence  difficulties. 
The  formulation  of  the  motion  analysis  problem  as  presented  in  this  report  is 
exact  and  automatically  casts  the  model  in  a linear  form  with  Gaussian  meas- 
urement perturbations.  Thus,  linear  estimation  techniques  apply,  and  the 
problem  of  algorithm  stability  normally  encountered  in  nonlinear  estimation  is 
avoided  and  improved  system  performance  is  achieved. 

The  modeling  and  theory  presented  provide  a solution  to  the  general 
motion  analysis  problem.  However,  this  report  is  devoted  almost  exclusively 
to  the  important  special  case  where  all  bearing  measurements  are  generated  by 
a single  moving  observer  as  illustrated  in  figure  2.  The  dynamic  model  that 
results  when  the  observer  is  also  traveling  at  a constant  velocity  represents  an 
unobservable  process.  5 ■ ^ In  a typical  situation,  observer  motion  is  often  re- 
stricted to  straight-line  segments  and  each  segment  is  termed  a "leg"  or 
"phase"  of  the  motion  analysis  problem.  It  is  well  known  that  only  a relative 
solution  (i.e. , estimates  normalized  by  the  slant  range  rs)  can.  be  obtained  on 
the  first  leg.  To  obtain  an  observable  process  for  the  purpose  of  estimating  the 
actual  value  of  the  target  state,  an  observer  maneuver  is  reqt  red.  Thus,  esti- 
mation of  target  state  from  bearing  observations  requires  a minimum  of  two 
legs  of  data.  This  observability  property  of  the  dynamic  process  results  in  a 
piecewise  convergence  of  the  state  estimate.  It  is  this  property  that  distin- 
guishes the  passive  bearings-only  problem  from  the  more  conventional  target 
trackers  and  motion  analysis  systems.'5-7 


Figure  2. 


Special  Case:  All  Bearing  Measurements  Generated  by  a 
Single  Moving  Observer 
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In  the  next  section,  a representation  of  the  target-observer  motion  process 
that  is  inherently  linear  in  form  is  developed;  this  departs  from  past  efforts 
where  linearization  about  the  nominal  or  estimated  target  state  is  required. 
Determination  of  the  target  position  and  velocity  is  then  directly  derived  from 
the  application  of  linear  estimation  theory.  Next,  the  properties  of  estimators 
for  "piecewise  observable"  systems  and  the  initialization  of  the  recursive  al- 
gorithm are  examined  in  detail.  The  practical  implementation  of  the  motion 
analysis  algorithm  is  then  discussed,  and  experimental  results  are  presented. 


MODELING  OF  THE  TARGET-OBSERVER  DYNAMICS 


A plot  of  the  target  and  observer  motion  was  shown  in  figure  2.  The 


discrete-time  equation  for 

by 

the  target  state  assuming  constant  velocity  is  given 

— 
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where  tg  is  the  time  increment  between  data  samples,  and  4>(k  + l,k)  is  the 
discrete  transition  matrix.  The  observer  motion  is  unrestricted  and  given  by 


rOx(k+C 

1 0| 

I rOx(k) 
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Equation  (2)  can  be  recast  as 


To 


xQ(k+l)  = 4>(k+l,  k)xQ(k)  + 


where 


aVOx(kl 

AVOy<k> 


(2) 


(3) 


AVo(k)  V0(k+1)  - VqOd  (4) 

is  the  incremental  change  in  the  observer  velocity.  The  bearing  to  the  target 
is  defined  bv  the  relation 


rTx(k)  - r0x(k)l/|rTv(k) 


rOv<k>l, 


(5) 


:i 


tan  fi  (k) 
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or,  equivalently, 

frTx(k)  " r0x(k>l  cos  $<k)  “ frTy(k)  - rQ  (k)]  sin  = °*  (6> 

The  bearing  measurements  viewed  by  the  observer  are  noise  corrupted;  that 
is, 

0m(k)  = ft(k)  +^(k),  (7) 

where  V (k)?  is  assumed  to  be  a purely  random  sequence  with  zero  mean  and 
variance  op  (k).  When  the  measured  bearing  is  used  in  place  of  the  true  bearing 
in  equation  (6),  this  relationship  becomes 

[rTx(k)  ~ rOx(k)|  cos  ^m(k)  " frTy(k)  _ rOy(k)I  sin  ^m^ 

{(i'Tx(k)  - r0x(k)]  cos  ft(k)  - [rTv(k)  - i'Qy(k)l  sin  ft  (k)  }cos  v (k) 

- < lrTx(k>  - r0x(k>l  sin  ^ *k)  + frTv^k>  " rOv(k)'  COS  ^ (k)  ^ Sin  V (k) 

= - rg(k)  sin  t’(k),  (8) 

where  the  last  step  results  from  substituting  equation  (6)  and  identifying  the 
quantity 

rg(k>  trTx(k)  " r0x(k'l  sin  ft  (k)  + frTv<k)  - r(^  (k)]  cos  ft  (k)  (9) 

as  the  slant  range  between  the  target  and  the  observer.  Rearranging  equation 
(8)  to  isolate  the  unknown  quantities  yields 

r0x(k)  cos  0m(k>  - ’’oyOO  sin  £|n<k> 

r^(k)  cos  £|ri(k)  - rT  (k)  sin  /?m(k)  + rg(k)  sin  ( k i , (10) 

where  the  left-hand  side  is  available  as  a measurement.  Defining 

H(k)  = H[0m(k)]  = [cos  £m(k)  -sin/?m(k)  0 0]  (11) 

and 

z(k)  = H(k)xQ(k)  (12) 

as  the  measurement  matrix  and  measurement  equation,  respectively,  equation 
(10)  assumes  the  compact  form 

/.(k)  H(k)xT(k)  + T)  (k),  (13) 


I 


is  only  positive  semi-definite  and  does  not  achieve  full  rank  until  an  observer 
maneuver  occurs  (see  appendix  A). 
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ESTIMATION  OF  TARGET  STATES 


The  bearings-only  motion  analysis  problem  is  stated  as  follows:  Given 
the  vector-valued  random  process  x(k)  and  observed  random  variables  Z(k) 

[z(k)  z(k-l). . . zfl)]',  find  an  estimate  x(k|k)  that  minimizes  the  expected  loss 
E{L(«r)!  = E {L[ ( f x (k)  - x(k)|  j)>.  Under  assumptions  that  are  consistent  with 
practical  constraints  and  desired  behavior,  the  optimal  estimate  is  given  by  the 
conditional  expectation8 

x(k|k)  E[x(k)|  Z (k)  1 . (18) 

If  the  random  processes  are  Gaussian.,  the  conditional  expectation  is  identical 
with  the  orthogonal  projection  of  x(k)  on  the  measurement  space  . ! » 10  How- 
ever, if  L(e)  ancj  the  estimates  are  restricted  to  be  a linear  combination 

of  the  measurement  Z(k),  then  the  optimal  estimate  is  also  the  orthogonal  pro- 
jection. Thus,  linear  estimation  is  bettered  by  nonlinear  estimation  only  if  the 
random  processes  are  non-Gaussian  and  then  only  if  third-order  and  higher 
probability  distributions  are  considered.  s For  the  bearings-only  motion  anal- 
ysis problem,  iz(k)  is  generally  assumed  Gaussian,  and  the  near  equality  of 
equation  (15)  holds  for  typical  noise  environments;  thus,  linear  estimation 
yields  the  optimal  estimate  for  essentially  all  performance  criteria.  A unified 
treatment  on  the  theory  of  linear  estimation  is  now  presented. 

LINEAR  ESTIMATION 

Given  the  dynamic  system 

x (k+ 1 ) <£(k+l,  k)x(k)  (19) 

and  measurement  equation 

z(k)  - H(k)x(k)  + T)  (k),  (20) 

the  vector  Z(k)  of  equation  (i7)  can  be  written  as 

Z(k)  = A(k)x(k)  + N(k),  (21) 

where 

H(k) 

H(k-l)[4>(k-l,k)] 

A(k)  - | . (22) 

H(l)(#(l,k)| 


6 


and 


N(k)=  [T)(k)T)(k-l)...TMl)]\  (23) 

The  penalty  or  loss  assigned  to  state  error  defined  by 

L = E { 1 1 x(k)  - x(k I k) 1 12) 

= trE  { (x(k)  - x(kjk)]  (x(k)  - x(k|k]'} 

= trP(k)  (24) 

serves  as  a measure  of  the  quality  of  the  estimate,  where  trf  • ] denotes  the 
trace  of  the  matrix  and  P(k)  is  the  error  covariance  matrix.  This  loss  func- 
tion is  minimized  when  the  state  error  (x(k)  - x(k|k)]  is  orthogonal  to  the 
available  data;  that  is, 

E {[x(k)  - x(k|k)]  Z'(k)}  = (0)  null  matrix.  (25) 

The  state  estimate  assumes  the  form 

x(k|k)  = B(k)Z(k),  (26) 

where  the  linear  combination  of  the  measurements  defined  by  the  matrix  B(k) 
is  to  be  determined.  Substituting  equations  (21)  and  (26)  into  equation  (25) 
yields 

Px(k)A'(k)  - B(k)[A(k)Px(k)A’(k)  + R(k)l  = [0],  (27) 

where 

R(k)  = E[N(k)N'(k)] 

| °77  <k)  0 

I 0 cr^(k-l)  0 

. 0 


! 0 

I 

and  where 

Px(k)  E[x(k)x'(k)[  $(k,0)Px(0)$’(k,0)  (29) 


0 


(28) 
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is  the  covariance  matrix  of  x(k)  without  any  measurement,  and  Px(0)  is  the 
a priori  covariance  of  the  unknown  vector  or  state.  From  equation  (27),  the 
optimal  value  for  the  matrix  B(k)  is 


B(k)  P (k)A'(k)[R (k)  + A(k)P  (k)A’(k)] 

X X 


= [P~V)  + A'(k)R“1(k)A(k)]"1A'(k)R_1(k), 


(30) 


where  the  last  step  results  from  the  matrix  inversion  identity  of  equation  (B-l) 
presented  in  appendix  B.  Note  that,  if  P^fo)  = o~l  and  R = [0]— or,  equivalently. 


R is  finite  and  P (o) 
x'  ' 


Ijm  Oq  I,  then 


B(k)  A’(k)[A(k)A’(k)l 


-1 


(31) 


is  the  pseudoinverse  of  A(k)  yielding  a minimum  norm  solution  for  the  unknown 
x(k)  for  the  underdetermined  case  (fewer  equations  or  measurements  than  un- 
knowns). 11  When  the  number  of  measurements  exceeds  the  number  of  un- 
knowns, the  inverse  in  equation  (31)  no  longer  exists  and  one  is  forced  to  edit 
the  number  of  measurements  to  achieve  a compatible  system.  If  noise  is 
present  (R  ^ (0]),  the  identity  in  equation  (30)  holds  and  all  the  measurements 
may  be  employed  in  extracting  the  estimates.  Note  that,  if  P”1  [0],  implying 

no  a priori  information,  the  relation 


B(k)  [ A ’ (k ) R 1 (k)A(k)]-1  A’(k)R_1(k) 


(32) 


results,  which  yields  the  weighted  least-squares  solution  for  the  overdetermined 
case  (more  measurements  than  unknowns).  10» 11  The  transition  between  cases 
is  given  by  the  general  form  of  equation  (30),  which  based  on  equations  (21)  and 
(20),  represents  a "generalized" inverse  of  the  matrix  A(k).  (The  relationship 
of  B(k)  to  the  standard  generalized  inverse  or  pseudoinverse,  and  sequential 
estimation  via  the  Kalman  algorithm  is  treated  in  the  next  section. ) 


With  the  optimal  estimator  defined  bv  equations  (20)  and  (30),  the  error 
covariance  matrix  in  equation  (24)  becomes 


P(k|  k)  P^(k)  - B(k)A(k)P^(k) 


Px(k)  - Px(k)A’(k)[R(k)  + A(k)Px(k)A,(k)f1A(k)Px(k), 


(33) 


where  the  first  step  utilizes  the  fact  that  the  optimal  value  of  x is  orthogonal 
to  the  state  error  and  K { [x-x|  (x-x)'}  E { (x-x)x' } . By  use  of  the  matrix 
inversion  lemma  equation  (B-2)— appendix  B,  the  optimal  error  covariance 
matrix  becomes 


P(k|k)  (P~ 1 (k ) 


A'(k)R  1 (k)  A (k ) |” 


(34) 
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Combining  equations  (34),  (30),  and  (26)  yields  as  the  optimal  estimate 

x(k|k)  = P(k)A'(k)R_1(k)Z(k).  (35) 

The  updating  of  the  state  estimate  x(k|k),  and  its  measure  of  quality 
P(k|k),  as  new  measurements  become  available  is  efficiently  accomplished  by 
a recursive  algorithm.  Recasting  equation  (34)  in  the  form 

P_1(k|k)  = P^(k)  + A'  (k)R~ 1 (k)A(k),  (36) 

utilizing  equations  (22),  (28),  and  (29),  and  noting  that 

IH(k+l) 

A(k+1)  = 

|A(k)<J>(k,k+l)_ 

permits  the  updating  of  the  error  covariance  matrix  to  be  expressed  in  the  form 
P_1(k+1  |k+l)  = [<l>(k+l,  k)P^(k)4>'(k+l,  k)]-1 

+ [$' 1 (k+1  ,k)]  'A'fl^R^fk)  A(k)$“ 1 (k+1 , k) 

+ -3—^ H'(k+l)H(k+l) 

cr“  (k+1) 

= [ $ (k+1 , k))_  1 [P~ 1 (k)  + A'(k)R_1(k)A(k)]<I>'1(k+l,k) 

+^ H’(k+l)H(k+l).  (37) 

CT^(k+l) 

Substituting  equation  (34)  into  equation  (37)  and  invoking  the  matrix  inversion 
lemma  (equation  (15-2 ) ) yields  the  computational  convenient  form 

P (k+1 1 k+1)  P(k+1  |k)  - P(k+1 1 k>H’ (k+1) 

• [H(k+l)P(k+l|  k)H’(k+l)  +a^(k+l)l_1 

• H(k+l)P(k+l|k),  (38) 

where 

P(k+l[k)  4>(k+l , k)P(k |k)  $ (k+1,  k).  (39) 
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Equation  (39)  represents  the  predicted  value  of  the  covariance  matrix,  and 
equation  (38),  the  corrected  or  updated  value.  Note  that  by  direct  expansion 
of  equation  (3(5)  it  can  be  shown  that 


P“Ylk)-  {3>(k,0)Px(0)  [^’(k.O)]}'1 
k 1 

+ L —±—  {[H(i)*(i,k)]'}  <[H(i)$(i,k)]}.  (40) 

i 1 an(i) 

Observability  assumes  that  the  second  term  of  equation  (40)  is  positive  definite. 
However,  as  noted  earlier  for  the  single  moving  observer  bearings-only  motion 
analysis  problem,  the  observability  of  the  matrix  depends  upon  whether  or  not 
an  observer  maneuver  has  occurred.  This  behavior  and  the  problem  of  initial- 
izing the  covariance  matrix  are  detailed  in  the  next  section. 


From  equation  (35)  the  updated  state  estimate  following  a new  measure- 
ment is  given  by 

x(k+l|k+l)  P(k+1  [k+l)A'(k+l)R_1(k+l)Z(k+l) 


P(k+1)  Jk+l)[H'(k+l)  j <£'(k-l,  k)A'(k)] 


2 1 

cr“  (k+1)' 
rj  ' i 


0 


' R(k) 


Z(k)  _ 


(41) 


After  some  manipulation,  the  substitution  of  equation  (38)  into  equation  (41) 
yields 


x(k+l J k+1 ) x(k+l|k) +K(k+l)(z(k+l)  - H(k+l)x(k+l|k)],  (42) 

where 

x(k+l|k)  = $(k+l,k)x(k|k)  (43) 


is  the  predicted  target  state,  and  the  correction  based  on  the  latest  measure- 
ment is  applied  with  gain 

K(k+1)  P(k+1  |k)H'(k+l)[H(k+l)P(k+l  |k)H'(k+l)+  cr^(k+l)|-1.  (44) 

Equations  (38),  (39),  (42),  (43),  and  (44)  are  the  Kalman- Buev  filter  equations. 
7'his  link  between  regressive  least  squares  and  Kalman  filtering  for  observable 
systems  was  established  earlier  by  Bryson  and  Ho1 and  Pagin'  *. 
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PROPERTIES  OF  LINEAR  ESTIMATORS  AND  INITIALIZATION 
OF  THE  RECURSIVE  ALGORITHM 

Given  the  set  of  simultaneous  linear  equations 

Z(k)  = Aik)x(k)  + N(k),  (21) 

where  without  loss  of  generality  it  is  assumed  that  the  data  are  normalized  such 
that  E[NN']  = a2I>  the  general  solution  is  classically  given  by15 

xLS(k)  " A (k)Z(k)  + [I  - A#(k)A(k)]W.  (45) 

The  matrix  A"(k)  is  the  pseudoinverse  of  the  matrix  A(k),  and  W is  an  arbitrary 
vector.  15  The  pseudoinverse  of  A(k)  is  defined  by 

A#(k)  =[A’(k)A(k)]“1A'(k)  (46) 

when  A'(k)A(k)  is  of  full  rank,  and  by 

A#(k)  [A'(k)A(k)]#A’(k)  (47) 

when  the  system  is  unobservable.  Since  A'(k)A(k)  is  symmetric,  an 
orthogonal  transformation  T exists  such  that 

A'(k)A(k)  = TAT',  (48) 

where 


0 

0 

• 

0 

0 

*2 

0 

• 

0 

0 

0 

• 

• 

1 • 

• 

A 

n~qA 

• 

• 

• 

• 

• 

0 

• 

, * 

* 

• * 

0 

0 

• 

• 

0 

diag[Aj)  , 

101 

1 

1 

(49) 

[01 

[01 
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where  n is  the  number  of  states  and  q ^ is  the  degeneracy  of  A'(k)A(k).  Thus, 
the  pseudoinverse  of  a symmetric  matrrx  is  defined  by 

[A'(k)A(k)l*  = T A T' 

r~ 

diag 

I 

= T 

[01 

The  value  of  x given  by  equation  (45)  represents  the  best  estimate,  in  a 
weighted  least-squares  sense,  to  equation  (21).  This  solution  is  unique  when 
A'(k)A(k)  is  nonsingular  since,  by  equation  (46),  the  second  term  in  equation 
(45)  vanishes.  When  the  process  is  unobservable,  that  is 

k 

det[A'(k)A(k)]  = det  £ [H(i)  $ (i,k)]’[H(i)  *(i,k)]  = 0,  (51) 

i=l 

a (q^)  parameter  family  of  solutions  for  x(k)  exists.  However,  even  in  this 
case*,  particular  types  of  solutions  can  be  singled  out  by  imposing  additional  con- 
straints. A common  selection,  from  the  infinity  of  solutions,  is  that  of  the 
minimum  norm.  Since  the  two  vectors  comprising  the  right-hand  side  of  equa- 
tion (45)  are  orthogonal,  10,15  the  minimum  norm  solution  occurs  when  W = [0], 
and  is  given  by 

x(k)  = A*  (k)Z(k).  (52) 

In  the  motion  analysis  problem  for  the  relative  target  state  x = x,p  - x , 
this  solution  (i.e. , equation  (52))  yields  the  null  vector  during  the  entire  first 
phase  (leg).  However,  under  these  conditions,  a range-normalized  solution  is 
often  sought.  From  the  foregoing  discussions,  it  is  clear  that  such  a relative 
motion  solution  must  be  derived  from  the  second  term  in  equation  (45). 

To  compute  Xj  ,_,(k)  sequentially,  one  employs 

xLg(k)  PLS(k)A'(k)Z(k),  <53> 

where,  in  this  case, 

PLg(k)  = (A’(k)A(k)]*.  <54> 

As  with  the  processing  of  equation  (34)  in  the  preceding  section,  use  of  the 
matrix  inversion  lemma  is  the  key  to  the  sequential  updating  of  x(k  |k)  as  new 
measurements  become  available.  However,  until  the  number  of  measurements 
equals  the  number  of  states  and  during  the  entire  first  leg  of  the  single  moving 
observer  problem  the  process  is  not  observable  and  A'(k)A(k)  is  singular. 
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Under  these  conditions,  algorithms  capable  of  generating  the  exact  pseudoin- 
verses are  given  by  Cline.  J-6  However,  in  this  report,  consideration  is  given 
to  the  use  of  a positive  matrix  P (k| 0)  in  computing  the  inverse,  that  is, 
approximate  PLg(k)  in  equation  (04)  by 


P(k|k)  = PLg(k)  = [Px(k|0)_1  + A'(k)A(k)]  1 

<!>  (k,  0)  [P-1(0  | 0 ) + <f>’  (k.O)A'(k)  A(k)$(k,0)]_V(k,0),  (55) 

where  P(0|0)  is  large.  Note  that  this  form  for  P(k|k)  is  directly  achieved  in 
the  development  presented  in  the  preceding  section  and  results  in  the  standard 
Kalman  computation,  with  P(0|0)  representing  the  initial  error  covariance 
matrix.  I he  effect  of  P(0|0)  on  the  computation  of  P(k)  for  the  specific  choice 

P(0|0)=crjjl  (56) 

is  now  examined. 


Substituting  equation  (56)  into  equation  (55)  yields 


P(k|k)  = <!>  (k,  0)  [ I + A^(k)A1(k)l  <£’  <k,  0),  (57) 

CT0 

where  A^(k)  = A(k)$(k,  0).  Since  A'  (k)A  (k)  is  symmetric,  an  orthogonal 
transformation  exists  such  that 

PLS1(k)  = Ai<k>Ai<k>  = Ti^iTr  (58) 

where 


Substituting  equation  (58)  into  equation  (57)  yields 


P(kjk)  = <*>(k,0)  T [±-  I +0,]  1T’®'(k,0) 

CTo  1 
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= * (k, 0)  T 


1 

n 

diag[— 9 + at.] 
_ 

1 

1 

1 

1 

[0] 

[0] 

1 

1 

f 

—o  1 

1 

1 

a0 

TNf’(k,0) 


$ (k,  0)  T 


drag  9 

+(T0 

, [oil 
« 
i 
1 

[01 

i 

1 9 

4>'  (k,  0). 


Separating  the  observable  and  unobservable  components  of  P(kjk)  permits 
equation  (59a)  to  be  recast  as 


r 

diag  - 


P(k'k)  = $(k,  0)  1 . 


or  “ <fxk,  0)  T 
0 1 


Tjd>'  (k,  0) 


T^4>’  (k,  0). 


Since- 

t 

r—  -i 

diag 

1 

= diagi 

1 1 1 

V“T 

- a0J 

1 2 
[Wi  °oaij 
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it  follows  that 

P(klk)  = $ (k,  0)$'  (k,  0)  [$  (k,  0)  (k,  0 + \ r*LS<k)]“1P1  g(k) 

°0 

+ cTq  * <k,  0)  (k,  0)  [ I - Pj'  g(k)PLS(k)], 

= PIS(k)  o P£sl*(k.0)*’(k,0)  + ^9  p:  g(k)]_1PLS(k) 
a0  * CT0 

+ o q <J>ik,0)4>’  (k,  0)  [1  - A (k)  A(k)l . (60) 


Substituting  equation  (60)  into  equation  (53)  shows  that  the  measurements  are 
indeed  processed  by  a matrix  that  approximates  A ” (k).  Specifically, 


P(k|k)A'(k) 


A (k)  = 


A (k)  -■ 


- PLS(k)[«(k,0)#,(k, 


0)  + 


V 


'0 


2 LS 
70 


(k)]_1A#(k), 


(61) 


where  use  of  the  identity  A AA'  = A'  has  been  made.  Thus,  the  introduction  of 
an  appropriate  nonsingular  P(0|0)  (specifically,  P(0|0)  as  detailed  above) 

permits  the  standard  recursive  algorithm  to  be  employed  in  the  calculation  of 
the  pseudoinverse  matrix  A'* (k).  The  second  term  in  equation  (60)  has  no  effect 
on  either  the  estimates  or  the  Kalman  gain  since 


P(k|k)H'  (k)  = $(k,0)*'(k,0)[$(k,0)*»(k,0)  + Pj  g(k)l  ^ g(k)H'(k). 

CT() 

I u 

The  effect  of  P(0|0)  on  the  estimates  clearly  depends  on  how  closely  A"  approx - 
imate§,A^.  Asa  cautionary  note,  it  should  be  pointed  out  that  in  selecting  values 
for  oq  for  initializing  the  computation  of  P(k),  the  choice  of  an  overly  large 
value  results  in  computational  difficulties  since  the  recursion  formulas  involve 
the  differencing  of  large  quantities. 

The  optimal  state  estimate  is  defined  by  equation  (42),  which  upon  re- 
grouping of  the  terms  becomes 

x(k  + l|k+l)  = [I  - K(k+l)ll(k+l)l<t><k+l, k)  x (k  k)  + K(k  + l)z(k+l).  (62) 
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The  general  solution  to  equation  (62)  is 

k 

x(k!k)  = 0 (k,  0)x(0|0)  <-  S 9 (k,j)K(j)z(j), 

i=l 


where 


0(k+l,k)  = [I  - K(k+l)H(k+l)]$>(k+l,k). 

It  is  shown  in  appendix  C that 
k 

r,  O (k,  j)K( j)z(j)  A (k)Z(k), 

j=l 

0 (k, 0)  [I  - A (k)A(k)14(k, 0), 
which  allows  equation  (63a)  to  be  recast  as 

x(k|k)  = A*(k)Z(k)  + [I  - A#(k)A(k)]x(k|0). 


(63a) 


(63b) 


(64) 


Thus,  when  the  process  is  unobservable,  the  minimum  norm  solution  is 
realized  by  initializing  the  recursive  algorithm  with  x(0|0)  = [0].  To  obtain 
other  solutions,  such  as  the  range-normalized  solution  during  the  first  leg  of 
the  single  moving  observer  problem,  requires  an  appropriate  selection  of 
x(0|0). 

The  relationship  of  the  estimate  x(k|k)  with  the  classical  least-squares 
solution  can  be  established  by  substituting  equation  (61)  into  equation  (64)  and 
utilizing  equation  (45),  where  W x(kjO),  Specifically, 


xlg(k) 


— , Pj  s(k)f$(k,0)$’(k,0) 
rr~ 


+ \ PLS(k)]  1rxLS(k)  - X(k|0)], 

a (65) 


where  for  observable  systems  |'n  Pj  ^(k)  (01  and  A (k)A(k)  = I.  lhus, 

the  estimate  of  the  Kalman-style  algorithm  asymptotically  yields  the  pseudo- 
inverse  solution. 
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APPLICATION  TO  THE  MOTION  ANALYSIS  PROBLEM 

This  section  is  concerned  with  various  state  representations  for  the 
motion  analysis  problem  and  the  resulting  estimation  algorithms.  For  the  tar- 
get state  xT,  the  results  of  the  preceding  section  apply  directly  and  the  esti- 
mate is  given  by 

xT(k+l[k+l)  = xT(k+l|k)  + K(k+l)[z(k+l)  - H(k+l)xT(k+l  [ k)]  , (66) 


where 


xT(k+l|k)  = $(k+l,k)xT(k]k). 

This  representation  is  preferred  for  the  general  problem  of  multiple  observers. 
However,  for  the  single  moving  observer  problem,  the  relative  target  state 

rx'k> 


Xr  = xt-  xO 


r (k) 
Vx(k) 

_v,«, 


is  often  of  interest  and  equation  (42)  reduces  to 

xR(k+l|k+l)  xR(k+l|k)  + K(k+l)[H(k+l)x0(k+l)  - h(k+l)xT(k+l|k)] 


xR(k+l  | k)  - K(k+l)H(k+l)xR(k+l[k), 


(67a) 


where  the  predicted  state  is  given  by 


xR(k+i|k)  = $(k+i,k)xR(k]k)  - 


A V q (k). 


(67b) 


The  estimated  target  velocity  vector  VT(k+l|k+l)  can  be  obtained  by  adding 
observer  velocity  components  to  Vx(k+l[k+l)  and  Vy(k+1  |k+l);  that  is. 


VTx(k+lJk+l) 

VTy(k+l  I k+1 ) 

The  relationship  between  the  various  estimates  is  illustrated  in  figures  3a 
through  3c,  where  the  covariance  and  gain  matrices  are  generated  via  equations 
(33)  and  (44),  respectively.  The  recursive  estimation  of  the  relative  target 
state  xR(k)  is  summarized  in  table  1. 


Vx(k+1 

k+1) 

\x(k+1) 

— 

A 

+ 

Vy  (k+1 

k+1) 

_V0y(k+n 
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Table  1.  Summary  of  the  Recursive  Estimation  of  the 
Relative  Target  State  xR(k) 


State  Vector: 


xR(k)=(rx(k)  ry(k)  Vx(k)  Vy(k)]' 
Initialization:* 


xR(0|0)  = [re(0)  sin)3m(0)  ro(0)  cos  0m(O)  0 0] 


m 


m 


P(0|0)=(TqI,  where  CT^  = constant 
Predicted  State: 

1)  0 I 


xR(k+l[k)  = <l>(k+l,k)xR(k[k)  - 


Predicted  Covariance  Matrix: 


0 0 0 


AVQ(k) 


P(k+1 1 k)  = $(k+l,k)P(k|k)$’(k+i,k) 
Measurement  Matrix: 


H(k+1)  = [cos  j3m(k+l)  -sin  )8  (k+1)  0 0] 

Gain  Matrix  Computation: 

K(k+1)  = P(k+1 1 k)H'(k+l)  [H(k+l)P(k+l  [ k)H'(k+l ) + r“(k+l)a2  (k+1)]"1 
Updated  State: 


xR(k+l  I k+1)  = xR(k+l  | k)  - K (k+1  )H (k+1  )xR (k+1  ] k) 

Updated  Covariance  Matrix: 

P(k+1  |k+l)  = [I  - K(k+l)H(k+l)lP(k+l|k) 


Target  Parameters: 


r8(k+l|k+l)  = (?x(k+l]k+l)  + ry  (k+1  [k+1)]1  2 
C-j-tk+l  [k+1)  - tan-1  [VTx(k+l  J k+1)/ V-pyfk+l  | k+1 11  “ 
ST(k+lJk+l)  = [V2  (k+ijk+l)  +y2  (k+i|k+l)I1/2 


Tv' 

A' 


£(k+l  J k+1)  tan"1  [rx(k+l  Jk+l)/r'(k+ 1 1 k+1 


♦Initialization  is  discussed  in  the  next  section. 


(1-1) 
<1—2 ) j 

(1-3) 

(1-4)| 

(1—5 ) j 

(1-6)  ^ 

i 

(1-7) I 

(1-8) 

(1-9) 

(1-10) 

(1-11) 

(1-12) 

(1-13) 
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IMPLEMENTATION  OE  THE  SINGLE  OBSERVER  TMA  ALGORITHM 


The  computation  of  the  gain  and  covariance  matrix  in  table  1 requires 
knowledge  of  the  slant  range  rs(k),  which  is  initially  unknown  to  the  observer. 
However,  from  equations  (52)  and  (53)  it  should  be  noted  that  lack  of  knowledge 
of  the  noise  variance 


Vk)~  rs(k)CTi;<k> 

in  the  computation  of  P(k)  does  not  preclude  achieving  an  estimate.  For  ex- 
ample, if  the  noise  variance  is  assumed  stationary,  representing  equal  weight 
on  all  data,  a least-squares  estimate  results.  Since  the  role  of  the  noise  var- 
iance in  the  computations  is  to  assign  a weight  to  each  measurement,  a near- 
optimal  implementation  can  be  realized  when  the  range  is  slowly  varying  by 
assigning  a constant  value  for  rs(k).  Introducing  the  normalized  covariance 
matrix  Pj^(k|k)  such  that 

P(kjk)  rg(k)PN(k|k)  (69) 

permits  equations  (1-5)  and  (1-9)  to  assume  the  normalized  form 


P^(k+1 1 k+1) 


rfo> 

r“(k+l ) 

s 


PN(k+l|k> 


PN(k+l|k)H'(k+l)H(k+l)PN(k+l[k) 


H(k+l)PN(k+l]k)H’(k+l)+a£(k+l) 


PN(k+l|k) 


PN(k+llk)H’(k+l)H(k+l)PN(k+lJk) 

H(k+l)PN(k+l|k)H'(k+l)+a£(k+l) 


(70) 


where 


PN(k+l|k)  4>  (k+1,  k)PN(k|  k)  <I>  (k+1,  k), 


and  where  a~  is  the  variance  of  the  noise  perturbing  the  bearing  measurements, 
and  the  last  step  in  equation  (70)  is  approximately  satisfied  for  slowly  varying 
rs(k).  Similarly,  for  the  gain  matrix,  substitution  of  equation  (69)  into  equation 
(1-7)  yields 


P (k+l]k)H'(k+l) 

K(k+1)  ■ l- ; o (71) 

H(k+l)P  (k+l|k)H'(k+l)+a^(k‘l) 
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Thus,  the  computation  of  both  P(k+i[k)  and  K(k+1)  is  independent  of  target 
range  when  P^(0|0)  is  chosen  independent  of  rg(0). 

In  the  initialization  of  the  state  vector,  it  is  seen  from  equation  (65)  that 
imposing  a minimum  norm  solution  on  the  relative  target  state  xp  xj  - Xq 
(i.e.  , setting  xp(()j  0)  [0 ) ) results  in  the  trivial  solution  xp(k|k)  (0|  for  all 
k on  the  first  leg.  To  generate  a relative  motion  (range-normalized)  solution, 
it  is  observed  from  equation  (1-8)  that  the  gain  matrix  operates  on  the 
"estimated'’  cross- range  residual;  that  is, 

H(k)x  (k+l|k)  r^(k-l  | k)  cos[ ft  (k- 1)  + v (k+1)]  - (k+ 1 J kt  sin[£(k+l)  • F(k«T)] 

-r(k+l]  k)  sin[  ft  (k+1 1 - ft  (k+1  j k)  + v (k- 1)] 

-r(k+l  j k ) sinf  Aft  (k+1)  + P(k-1)].  (72) 


Thus,  if  in  the  initial  state  estimate  the  relative  velocity  is  zero,  the  solution 
for  xp(k)  is  proportional  to  (or  scaled  bv)  the  initial  slant  range  estimate. 
This  is  a desirable  characteristic  for  the  estimator,  and  the  initialization  of 
the  state  estimate  given  by 

r (0)  sin  ft  (On 
s m 


xR(0|0) 


? (0)  cos  ft  (0) 
s m 


is  used  in  this  report.  On  the  second  leg,  the  svstem  becomes  observable; 
thus,  the  second  term  in  equation  (45)  fades  and  the  significance  of  the  initial 
value  for  xp  diminishes.  In  the  initialization  of  the  covariance  matrix,  exper- 
imental studies  have  shown  that  the  choice 


P(0 1 0)  ml, 


2 2 

where  (0 ) , produces  good  results  in  realizing  the  desired  estimator 

performance.  From  equation  (69)  i)  is  seen  that  setting  r~(0)  is  equiv 

lent  to  initializing 


equiva- 


PN(0[0) 


This  initialization  can  be  shown  to  generate  weights  on  new  data  in  a manner  that 
takes  into  account  the  size  of  the  bearing  change  between  samples  in  relation  to 
the  bearing  error  due  to  measurement  noise.  ' ' The  normalized  algorithm,  in- 
cluding initialization,  is  summarized  in  table  2. 
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Table  2.  Summary  of  the  Normalized  Recursive  Algorithm 
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PROPERTIES  OF  THE  TMA  SOLUTION 

The  convergence  of  the  recursive  solution  (Kalman)  to  the  least-squares 
estimate  as  the  number  of  measurements  increases  — see  equation  (65)  — 
permits  the  behavior  of  the  linear  TMA  estimator  to  be  discussed  in  terms  of 
the  idealized  estimate  of  equation  (45).  Note  that  the  measurement  matrix  A(k), 
defined  by  equation  (22),  is  a function  of  the  measured  bearing  and  hence  the 
measurement  noise  p(k).  It  is  shown  in  appendix  A that  under  noise-free 
conditions  A'(k)Aik)  is  singular  over  the  entire  first  leg.  Since  for  the  TMA 
problem  the  minimum  norm  solution  x(k)  = A“Z(k)  is  the  null  vector,  the  first 
leg  solution  prior  to  an  observer  maneuver  must  be  the  second  term  of  equation 
(45).  This  represents  the  range-normalized  solution. 

In  the  experiments  described  in  the  next  section  this  solution  is  indeed 
observed  under  low-noise  conditions.  However,  when  the  level  of  the  measure- 
ment noise  is  substantial,  A'(k)A(k)  becomes  nonsingular  even  on  wit  Hrst  leg; 
and  hence  the  factor  [I  - A#(k)A(k)]  in  equation  (45)  vanishes  and  the  solution 
for  any  initial  conditions  collapses  toward  the  minimum  norm  solution.  This 
behavior  appears  to  be  abated  in  the  low  -noise  case  bv  preprocessing  of  the 
data  and  approximate  computation  of  Aw(k).  To  demonstrate  the  consistency  of 
the  first  leg  behavior  when  A'(k)A(k)  is  nonsingular,  i.  e. , to  show  that  the 
addition  of  noise  to  the  matrix  A(k)  does  not  alter  system  observability,  recall 
that  the  measurements  are  generated  from  the  measured  bearings  and  observer 
position  by  equation  (12);  that  is 

z(k)  = H(k)x  (k)  H(k)xT(k)  + r)(k). 

On  the  first  leg, 

Z(k)  - (z(k),  z(k-l) z(l)]'  A<k)x()(k),  (76) 

since  observer  velocity  is  constant.  Consequently,  the  estimate  is  given  by 
x(k|k)  = A#(k)Z(k)  + [I-A*(k)A(k)]x(kjO) 

= A"(k)A(k)xT(k)  + A"(k)N(k)  + (I  - A"(k)A(k)]x(k|0) 

A*(k)A(k)xQ(k)  [l  - A"(k)A(k)]x(ki0);  (77) 

or,  alternatively,  the  error  on  the  first  leg  is 

^ % 
xe(k|k)=  [I  - A#(k)A(k)][xT(k)  - xT(k  O)]  - A '(k)N(k) 
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= xT(k)  - A#(k)A(k)xQ(k)  - [I  - A#(k)A(k)]xT(k|0).  (78) 


When  the  level  of  noise  is  substantial,  A'(k)A(k)  becomes  nonsingular;  and 

xg(k|k)  = -A  (k)N(k)  = xT(k)  - x0(k).  (79) 

Thus,  prior  to  an  observer  maneuver,  it  is  not  possible  to  identify  the  complete 
target  state  and  the  TMA  estimate  tends  toward  the  minimum  norm  solution. 
Following  a maneuver,  equation  (69)  has  an  additional  term  and  the  estimate 
approaches  the  actual  target  state.  The  effect  on  the  solution  generated  by  the 
algorithm  of  table  2 is  given  by  equation  (65)  and  can  be  shown  to  produce  a bias 
that  fades  as  the  number  of  measurements  increases.  Substituting  equation  (21) 
into  equation  (45)  yields 

xT(k|k)  = xT(k)  + A (k)N(k),  (80) 

which  reveals  an  additional  source  of  bias.  Since  both  A (k)  and  N(k)  are 
functions  of  the  measurement  noise  p(k),  the  expectation  of  the  second  term  is 
non-zero.  This  bias  is  structural  and  not  inherent  in  nature;  it  is  a common 
occurrence  in  parameter  estimation  problems  in  control  theory,  and  in  the 
linear  prediction  techniques  of  signal  processing.  19 
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EXPERIMENTAL  RESULTS  I 


Representative  target-observer  geometries  were  simulated  on  the  AN/ 
UYK-7  computer  in  the  Systems  Analysis  Laboratory  at  the  Naval  Underwater 
Systems  Center.  White  Gaussian  noise  with  zero  mean  and  0.5  standard  devia- 
tion was  added  to  the  data  generated  by  the  dynamic  model  as  shown  in  figure  4. 
Preprocessing  of  the  target  bearing  was  simulated  by  averaging  twenty  1 -second 
samples  to  provide  a bearing  measurement  for  a single  (moving)  observer  tar- 
get motion  analysis  algorithm. 


Figure  4.  Simulation  Diagram  of  Motion  Analysis  System 

The  vehicle  tracks  for  the  first  target-observer  geometry  considered  are 
shown  in  figure  5a.  Both  vehicles  are  traveling  at  a constant  speed  of  5.626 
meters/sccond,  and  the  observer  traverses  four  240-second  legs  involving 
course  changes  of  90°.  In  this  first  example,  instantaneous  or  point  maneuvers 
are  used. 

The  relative  motion  plots  of  the  solutions  obtained  on  the  first  leg  (shown 
in  figure  5b)  demonstrate  the  property  of  the  algorithm  (summarized  in  table  2) 
to  provide  solutions  scaled  by  the  initial  range  estimate  under  the  conditions  of 
low  measurement  noise  and  no  jitter  in  observer  velocity.  Figure  5c  shows 
the  true  and  estimated  target  tracks  for  various  initial  range  estimates.  Note 
that,  as  predicted,  the  dependency  of  the  solution  on  the  initial  range  estimate 
fades  rapidly  once  the  process  becomes  observable.  Figure  5d  through  5g  show 
the  behavior  of  the  components  of  the  estimated  target  state  £^(k|k)  as  a function 
of  elapsed  time.  It  should  be  noted  that,  on  the  second  leg  for  the  particular 
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geometry  considered,  the  vehicles  are  running  parallel  and  the  rate  of  change  of 
the  bearing  angle  is  zero.  Under  these  conditions,  except  for  the  transient 
immediately  following  the  maneuver,  no  further  improvement  in  the  velocity 
estimate  along  the  line  of  sight  is  possible  and  this  is  reflected  in  figures  5d 
and  5f. 

The  target  motion  parameters  often  of  interest  are  the  range,  course, 
and  speed.  Plots  of  the  range,  course,  and  speed  errors  for  the  first  example 
are  shown  in  figures  6a  through  6c  for  two  initializations  (r(0)  10,  000  meters 

and  r(0)  0,  which  is  the  minimum  norm  solution).  The  errors  via  the  appli- 

cation of  the  extended  Kalman  filter  (based  on  the  linearization  of  equation  (5) 
about  a nominal  trajectory,  see  appendix  D)  are  also  shown  for  purposes  of 
comparison.  Note  the  substantial  course  error  and  partial  divergence  of  the 
range  solution  that  occurs  on  the  second  leg  when  the  latter  method  is  employed. 
Fortunately,  in  this  example,  the  situation  is  corrected  following  the  next  ob- 
server maneuver. 

Two  other  encounters  are  shown  in  figures  7a  through  7d  and  figures  8a 
through  8d.  In  both  of  these  cases,  the  turning  rate  of  the  observer  is  con- 
strained to  3°/second.  This  constraint  represents  the  effect  of  finite  thrust  in 
space  applications  or  the  restricted  turning  rate  of  hydrodynamic  vehicles  in 
marine  applications.  In  the  second  target-observer  geometry  (figure  7a),  the 
extended  Kalman  filter  exhibited  premature  collapse  of  the  covariance  matrix 
resulting  in  solution  divergence.  Such  instability  of  the  extended  Kalman  filter 
solution  in  the  target  motion  analysis  (T\IA)  problem  was  first  observed  by 
Fagin  (see  Murphy-).  This  instability,  however,  is  fundamental  to  nonlinear 
modeling  and  is  aggravated  in  this  particular  application  bv  the  fact  that  the 
estimated  state  used  in  the  linearization  process  is  not  observable  during  the 
critical  initial  period  of  the  problem.  ^ The  modeling  presented  in  this  report 
precludes  this  difficulty.  Bias  errors  in  the  estimate  did  not  appear  significant 
for  the  particular  geometries  and  noise  environment  examined. 
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Figure  5d.  Actual  and  Estimated  Relative  X-Position  (rx)  as  a 
Function  of  Elapsed  Time  (Geometry  No.  ll 
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Figure  6b.  Target  Course  Error  for  Various  Initialization  Parameters 

(Geometry  No.  1) 
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O — NEW  ALGORITHM  (SEE  TABLE  2) 

iR(0)  = [f  ,(0)  r‘y(0)  0 0 

rs  (Oi  = 10,000  METERS 

EXTENOED  KALMAN  FILTER 
(SEE  APPENDIX  D) 

* R (0)  = [f  * (0)  fy(0)  -V0x(0)  -v0y 

fs  (0)  = 10,000  METERS 


COURSE  ERROR (DEGREES) 
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O — NEW  ALGORITHM (SEE  TABLE  2 ) 

*r(0)  = [f , (0)  r*y(0)  0 0 

fs  (0)  = 10,000  METERS 

— EXTENDED  KALMAN  FILTE  R ALGORITHM 
(SEE  APPENDIX  D)  f 
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Figure  7c.  Target  Course  Frror  for  Geometry  No.  2 
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Figure  7d.  Target  Speed  Error  for  Geometry  No.  2 
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SUMMARY  AND  CONCLUSIONS 


Bearings-only  motion  analysis  is  modeled  in  the  form  of  a linear  state 
estimation  problem  when  the  measurements  are  generated  from  the  bearing  ob- 
servations via  equation  (12).  Implementation  of  equation  (12)  requires  knowl- 
edge of  the  emitter's  bearing,  location  of  the  observer,  and  time  of  observation. 
For  the  special  case  of  a single  moving  observer,  emphasized  in  this  report, 
the  resulting  dynamic  process  is  unobservable  on  the  first  leg.  The  relationship 
between  recursive  least-squares  and  Kalman  filtering  when  the  dynamic  process 
is  unobservable  is  detailed.  An  estimation  algorithm  (table  2),  including  initial- 
ization, is  presented  that  exhibits  consistent  behavior  for  all  iarget-observer 
geometries  examined.  Minimum  norm  solutions  are  provided  on  the  first  leg  of 
the  problem;  and  convergence  to  the  complete  solution  is  achieved  following 
maneuvers  by  the  observer  that  result  in  the  dynamic  process  becoming  observ- 
able. Convergence  of  the  estimate  of  target  state  is  dependent  upon  the  velocity 
variation  undertaken  by  the  observing  vehicle.  No  attempt  was  made  in  this 
report  to  optimize  the  observer  maneuvers  to  improve  convergence.  This  is  an 
area  fo^  further  study. 

Since  the  gain  and  covariance  matrices  are  dependent  only  upon  the  meas- 
ured bearings,  changes  in  target  velocity  have  the  same  effect  on  these  quanti- 
ties as  observer  maneuvers.  Thus,  evasive  action  by  the  target  can  produce 
false  solutions  and  convergence  of  the  covariance  matrix.  Random  perturbation 
in  the  target  velocity  can  be  handled  by  the  addition  of  plant  noise  to  the  model 
of  equation  il).  When  the  heading  changes  undertaken  by  the  observer  are  suf- 
ficiently large,  the  angle  tracking  may  require  new  antenna  systems;  and  the 
problem  is  further  complicated  by  the  addition  of  bias  errors  introduced  be- 
tween the  data  legs.  The  more  general  three-dimensional  problem  using  addi- 
tional passive  sensor  inputs  is  currently  under  study. 
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APPENDIX  A 

OBSERVABILITY  OF  THE  NOISE-FREE  TARGET-OBSERVER  PROCESS 


To  provide  insight  to  the  behavior  of  the  bearings-only  motion  analysis 
problem,  the  observability  of  the  noise-free  ease  is  examined.  Observability 
for  time-varying  systems  requires  that  the  n x n matrix 


2 (H(i)$  (i,k)]'[H(i)$(i,k)]  = A'(k)A(k) 

i 1 


(A—  1) 


be  positive  definite.  Since  the  rank  of  the  product  of  matrices  cannot  exceed 
the  rank  of  individual  matrices,  that  is, 


Rank  (AB)  ^ Rank  (A)  or  Rank  (B), 
Rank  (A 'A)  = Rank  (A), 


(A-2) 
(A- 3) 


the  system  is  unobservable  if  the  rank  of  A(k)  is  less  than  the  number  of  states, 
which  is  four  in  the  two-dimensional  motion  analysis  problem.  From  equations 
(11,  (11),  and  (22)  we  have 

H(k) 

H(k-l)<J>  (k-l,k) 

H(k-2)<J>(k-2,k) 


A(k) 


H(l)<i>(l,k) 
cos  0 (k) 
cos  0 (k-1) 
cos  0 (k-2) 


-sin  0 (k) 
-sin  0 (k-1) 
-sin  0 (k-2) 


0 

-t  cos  0 (k-1) 


t sin  $(k-l) 

o s 

-2ts  cos  0 (k-2)  2tssin  0(k-2) 


cos  0 (1) 


-sin  0 (1) 


(k-1  )tscos/S(lt  (k-1  (t^sin  0(1) 


. (A-4) 


A-l 


W 
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If  the  bearing  rate  is  zero,  implying  /3  (k+1)  8(k),  then  the  first  two  columns 

and  the  last  two  columns  of  the  matrix  in  equation  (A-4)  are  linearly  dependent 
and  the  rank  of  A(k)  is  two.  Physically,  this  implies  that  it  is  not  possible  to 
determine  the  slant  range  of  the  target  and  the  component  of  the  relative  velocity 
Vr  = V-p  - Vo  along  the  line  of  sight.  For  the  case  where  the  single  moving 
observer  is  traveling  at  constant  velocity  and  d/3/dt  ^ 0,the  rank  of  A(k)  can 
be  shown  to  be  three.  To  determine  that  the  matrix  A(k)  does  not  achieve  full 
rank  prior  to  observer  maneuvers,  note  from  equations  (12)  and  (13)  that  in 
the  noise-free  case 


H(k)xR(k)  (0),  (A-5) 

where  the  behavior  of  the  relative  target  state  (constant  observer  velocity)  is 
defined  by 


xR(k+l ) <f>  (k+1,  k)xR(k).  (A-6) 

By  comparing  equations  (A-5)  and  (A-6>  with  equations  (19)  and  (20),  one  ob- 
tains an  equation  that  is  similar  to  equation  (21);  that  is, 

A(k)xR(k)  f0[,  (A- 7) 

or 

A'(k)A(k)xR(k>  1 0 1 . (A-8) 

It  is  clear  that  in  order  for  equation  (A-mi  to  be  satisfied,  either  the  trivial 
solution  XR(k)  (0[  is  unique  or'  \(k>  is  not  of  full  rank,  in  which  case  non- 
trivial solutions  \R(k)  t 0 exist.  However,  from  equation  (5)  one  has 


8 (t> 


tan-! 


Vxt+rx(0) 

Vvt+rv<°> 


(A-9) 


and  it  is  clear  that  a one-parameter  family  of  solutions,  scaled  by  rs(0), 
exists.  This  situation  is  illustrated  in  figure  A-l.  Consequently,  the  rank  of 
the  matrix  A'(k)A(kl  is  not  equal  to  four,  and  the  system  is  not  observable 
during  the  first  phase  (leg)  of  the  bearings-oniv  motion  analysis  problem.  For 
the  system  to  become  observable,  an  additional  linearly  independent  measure- 
ment is  required;  this  is  realized  by  a change  in  the  velocity  of  the  observer 
(i.e. , by  an  observer  maneuver). 


A- 2 
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APPENDIX  B 

MATRIX  INVERSION  IDENTITIES 

Given  the  n x n positive  definite  matrix  P,  the  m x n matrix  A,  and 
the  m x m positive  definite  matrix  R,  then 

PA'fR  + APAT1  [P-1  + A'R"1A]"1A'R_1,  (B-l) 

and 

[P-1  + A’R-1A]“1  P - PA'  [R  + APA']-1  AP.  (B-2) 

Proof: 

1.  Pre-  and  post-multiply  equation  (B-l)  by  the  nonsingular  matrices 
[p-1  + A'R-1A]  and  (R  + APA'],  respectively. 

2.  Substitute  equation  (B-l)  into  (B-2),  yielding, 

[P_1  + A'R_1A]_1  P - (P_1  + A»R_1A]_1A»R"1AP, 
and  pre-multiply  by  the  nonsingular  matrix  fP  1 + A’R  1 A ] . 


B-l/ B-2 
Reverse  Blank 
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APPENDIX  C 

RECURSIVE  COMPUTATION  OF  THE  MINIMUM  NORM  SOLUTION 
Given  the  state  transition  matrix 

0(k+l,k)  |I  - K(k+l)H(k+l)]<l>  (k+l,k)  (C-l) 

of  the  optimal  estimator,  the  optimal  gain  matrix  K(k),  and  the  history  of 
measurements  Z (k)  [z(k)  z(k-l). . . z(l)]',  then 

k "V 

L 0 (k,  j)K(j)z(j)  A (k)Z(k),  and 
j=  1 rv/ 

0(k,0)  = [I  - A#(k)A(k)]4>(k,0).  (C-2) 

Proof: 

1.  Rewrite  equation  (C-l)  in  the  form 

0(k+l,k)  - (I-  K(k+1  )H(k+l  >|P(k+l  J k > P~ 1 (k+1  ]k)  <i>  (k+1,  k).  (C-3) 


Utilizing  the  relations 


P(k+1 1 k+1) 

[I  - K(k+1  )H(k+l)JP(k+l  | k), 

(C— 4 ) 

(k+1 , k) 

$(k,k+l), 

( C— 5 > 

and  . 

P"  (k+1  J k) 

<J> 1 (k,  k+1  )P_  1 (k  I k)  (k,  k+1), 

(C-6) 

the  optimal  state  transition  matrix  can  be  expressed  as 

0 (k+1, k) 

P(k+1  |k+l)  <f>'(k,  k+l)P  1 (k  (k), 

(C-7) 

or 

0 (k, j) 

P(k|k)4>'(j.  k)P_1(jJj), 

(C-8) 

for  all  values  of  j.  Since 

K(j)  = P(j|j)H'(j)/ar , 

(C-9) 

C-l 


w 


n 
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the  product  0 (k,j)K(j)  is  given  by 

0(k,j)K(j)  - P(k!  k)  k)P_1(j  | j)P0 1 j>H’(j)/crJ 


P(kjk)4>'(j,  k)H'(j)  o 


J 


(C-10) 


Substituting  equation  (C-10)  into  the  left-hand  side  of  equation  ( C— 2 ) ;md  ex- 
panding yields 


£ 0 (k,  j)K(j)z(j)  P(k]k)  H -f  ^,(.k  k.tH'fk-1  i/(k  li 


j=l 


a, 


k-l 


$'(k-2,  k)H'(k-2)/.(k-2)  «5' 1 1 , k)ll'( l)z( 1 ) 


k-2 


= P(kjk)A'(k)Z(k). 


(C-ll) 


Substitution  of  equation  (61)  into  equation  (C-lli  completes  the  proof. 

2.  Setting  j = 0 in  equation  (C-8)  yields 

0(k,  0)  = P (k , k><5» ' (0,k)P-1(0  0).  (C-12) 

Since  P ^(0|0)  = <J>'(k,0)P^  (k  0)$(k,0),  it  follows  that 
0(k,  0)  = PfkjklP”1  (k  |0)<J>(k,  0) 

= Jp(k  kHP"1  (k!0)  + A'(k)A(k)]  - P(kk)A’(k)A(k/$(k,0). 

J <C-13> 

Substituting  equation  (55)  into  equation  (C-13)  and  recognizing  that  A ' = 
P(kjk)A'(k)  yields 

0fk,  0)  = fi  - A (k)A(k)]  <J>(k,  0).  (C-14) 


C-2 
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APPENDIX  D 

SUMMARY  OF  THE  APPLICATION  OF  THE  EXTENDED 
KALMAN  FILTER  TO  THE  MOTION  ANALYSIS  PROBLEM 


A summary  of  the  application  of  the  extended  Kalman  filter  to  the  bearings- 
only  motion  analysis  problem  is  presented  as  follows: 

State  Vector: 


x(k)  [r^ik)  r (k)  V^(k)  Vv(kl]' 


(D—  1) 


Initialization: 


x ( 0 ; 0 ) [?g(0)  sin  0m(O)  rs(0)  cos /3m(0)  -V0x<k)  -VOv(k'*'  (U-2' 


Prediction: 


P(0j0)  = 


0 0 


0 0 


o o of, 

V v 


x(k+l|k)  <$  (k+1, k)x(k I k) 


(D-3) 


P(k+l|k)  <f>  (k+l,k)P(k|k)<J>'(k+l,k) 
/3(k+l | k ) tan  1 [r^(k+l  |k)/r^(k+l  [ k)] 

r (k+1 1 k > [r~(k+l  |k)  + f“(k+l ] k) ] 1 


(D— 5 ) 
(D-6) 


(D— 7) 


Measurement  Matrix  Linearization: 


1 1 ( k + 1 ) 


5/?  Z8_ 


10 

VTy_ 


x(k+li  x(k+l  k) 


( D— 8 > 


where 


cos  #<k+l | k)/r  (k+1 ] k> 


x(k+l ) x(k+l [ k) 


(D-iH 


I)-  I 


p 
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-sin  0(k+l|k)/r  (k+l|k) 

x(k+l)  x(k+l|k> 


0 


(D-10) 


(D-ll) 


TT7—  o (D-12) 

Ty 

Gain  Computation: 

K(k+1)  P(k+l|k)H'(k+l)fH(k+l)P(k+l|k)H'(k+l)  + CT^f 1 (D-13) 


Updating  of  Estimate: 

x(k+l  ik+1)  = x(k+l|k)+ K[0m(k+1)  - 0(k+l|k)]  (D-14) 

P(k+1  |k+l)  [I  - K(k-d  )H(k+l)]P(k+l  | k) 

= P(k+1 1 k)  - P(k+1 1 k)H'(k+l)H(k+l)P(k+l  | k) 

• [H(k+l)P(k+l  ]k)H’(k+l)  + aj  (k+1)]-1  (D-15) 

Estimated  Target  Parameters: 

Target  Range  rg(k+l|k+l)  [r~(k+l|k+l)  + r“ (k+1  J k-»- 1 ) 1 1 “ (D-16i 

a . _i  a . a 

Target  Course  C,j,(k4l|k+1)  tan  [VTx(k+l  | k+l)/VT^(k+l  | k4 1)]  (O- 1 7> 

Target  Speed  = ST(k+l|k+l)  [V^x(k+l|  k+1)  + V^.fk+l  | k+ 1 )1 1 (D-18) 


D-2 
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